Question: “How can I crop raster objects to vector objects, and extract the summary of raster pixels?”

Objectives:

  • Crop a raster to the extent of a vector layer.
  • Extract values from a raster that correspond to a vector file overlay.

Keypoints:

  • Use the crop() function to crop a raster object.
  • Use the extract() function to extract pixels from a raster object that fall within a particular extent boundary.
  • Use the extent() function to define an extent.
library(dplyr)
library(sf)
library(tibble)
library(ggplot2)
library(raster)
library(rgdal)
library(ggplot2)
library(dplyr)
library(here)

Introduction

copied from the carpentry lesson Manipulating Raster Data).

We often work with spatial layers that have different spatial extents. The spatial extent of a shapefile or R spatial object represents the geographic “edge” or location that is the furthest north, south east and west. Thus is represents the overall geographic coverage of the spatial object.

Image Source: National Ecological Observatory Network (NEON) The graphic below illustrates the extent of several of the spatial layers that we have worked with in this workshop:

Image Source: DCC

Frequent use cases of cropping a raster file include reducing file size and creating maps. Sometimes we have a raster file that is much larger than our study area or area of interest. It is often more efficient to crop the raster to the extent of our study area to reduce file sizes as we process our data. Cropping a raster can also be useful when creating pretty maps so that the raster layer matches the extent of the desired vector layers.

Import the raster

Data available here.

DSM_TUD <- raster(here("data","tud-dsm.tif"))
DTM_TUD <- raster(here("data","tud-dtm.tif"))
CHM_TUD <- DSM_TUD - DTM_TUD

CHM_TUD_df <- as.data.frame(CHM_TUD, xy = TRUE)

oai_boundary_tudlib <- st_as_sfc(st_bbox(raster(here("data","tudlib-rgb.tif"))))

Crop a Raster Using Vector Extent

We can use the crop() function to crop a raster to the extent of another spatial object. To do this, we need to specify the raster to be cropped and the spatial object that will be used to crop the raster. R will use the extent of the spatial object as the cropping boundary.

To illustrate this, we will crop the Canopy Height Model (CHM) to only include the area of interest (AOI). Let’s start by plotting the full extent of the CHM data and overlay where the AOI falls within it. The boundaries of the AOI will be colored blue, and we use fill = NA to make the area transparent.

ggplot() +
  geom_raster(data = CHM_TUD_df, aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
  geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
  coord_sf()

Now that we have visualized the area of the CHM we want to subset, we can perform the cropping operation. We are going to use the crop() function from the raster package to create a new object with only the portion of the CHM data that falls within the boundaries of the AOI.

CHM_TUD_Cropped <- crop(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib))

Now we can plot the cropped CHM data, along with a boundary box showing the full CHM extent. However, remember, since this is raster data, we need to convert to a data frame in order to plot using ggplot. To get the boundary box from CHM, the st_bbox() will extract the 4 corners of the rectangle that encompass all the features contained in this object. The st_as_sfc() converts these 4 coordinates into a polygon that we can plot:

CHM_TUD_Cropped_df <- as.data.frame(CHM_TUD_Cropped, xy = TRUE)

ggplot() +
  geom_sf(data = st_as_sfc(st_bbox(CHM_TUD)), fill = "green",
          color = "green", alpha = .2) +  
  geom_raster(data = CHM_TUD_Cropped_df,
              aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  coord_sf()

The plot above shows that the full CHM extent (plotted in green) is much larger than the resulting cropped raster. Our new cropped CHM now has the same extent as the aoi_boundary_HARV object that was used as a crop extent (blue border below).

ggplot() +
  geom_raster(data = CHM_TUD_Cropped_df,
              aes(x = x, y = y, fill = layer)) + 
  geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  coord_sf()

We can look at the extent of all of our other objects for this field site.

st_bbox(CHM_TUD)
    xmin     ymin     xmax     ymax 
 83569.5 445251.5  87180.0 447180.0 
st_bbox(CHM_TUD_Cropped)
    xmin     ymin     xmax     ymax 
 85272.0 446295.0  85661.5 446694.0 
st_bbox(oai_boundary_tudlib)
     xmin      ymin      xmax      ymax 
 85272.00 446295.20  85661.28 446694.24 
# plot_locations_HARV <-
#   read.csv("data/NEON-DS-Site-Layout-Files/HARV/HARV_PlotLocations.csv")
# point_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
# utm18nCRS <- st_crs(point_HARV)
# plot_locations_sp_HARV <- st_as_sf(plot_locations_HARV, coords = c("easting", "northing"), crs = utm18nCRS)
# 
# st_bbox(plot_locations_sp_HARV)

leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) %>% 
  filter(leisure %in% c("playground", "picnic_table"))
Reading layer `delft-leisure' from data source 
  `/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD New
st_bbox(leisure_locations_selection)
     xmin      ymin      xmax      ymax 
 81863.21 442792.82  86719.87 449007.92 

Our plot location extent is not the largest but is larger than the AOI Boundary. It would be nice to see our vegetation plot locations plotted on top of the Canopy Height Model information.

Challenge :

Solution

CHM_plots_TUDcrop <- crop(x = CHM_TUD, y = leisure_locations_selection)

CHM_plots_TUDcrop_df <- as.data.frame(CHM_plots_TUDcrop, xy = TRUE)

ggplot() + 
  geom_raster(data = CHM_plots_TUDcrop_df, aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  geom_sf(data = leisure_locations_selection) + 
  coord_sf()

In the plot above, created in the challenge, all the vegetation plot locations (black dots) appear on the Canopy Height Model raster layer except for one. One is situated on the blank space to the left of the map. Why?

A modification of the first figure in this episode is below, showing the relative extents of all the spatial objects. Notice that the extent for our vegetation plot layer (black) extends further west than the extent of our CHM raster (bright green). The crop() function will make a raster extent smaller, it will not expand the extent in areas where there are no data. Thus, the extent of our vegetation plot layer will still extend further west than the extent of our (cropped) raster data (dark green).

Image Source: DCC # Define an extent

So far, we have used a shapefile to crop the extent of a raster dataset. Alternatively, we can also the extent() function to define an extent to be used as a cropping boundary. This creates a new object of class extent. Here we will provide the extent() function our xmin, xmax, ymin, and ymax (in that order).

# extent(CHM_TUD_Cropped_df)
new_extent <- extent(85272.25, 85661.25, 446295.2, 446693.8)
class(new_extent)
[1] "Extent"
attr(,"package")
[1] "raster"

TIP: The extent can be created from a numeric vector (as shown above), a matrix, or a list. For more details see the extent() function help file (?raster::extent).

Once we have defined our new extent, we can use the crop() function to crop our raster to this extent object.

CHM_TUD_manual_cropped <- crop(x = CHM_TUD, y = new_extent)

To plot this data using ggplot() we need to convert it to a dataframe.

CHM_TUD_manual_cropped_df <- as.data.frame(CHM_TUD_manual_cropped, xy = TRUE)

Now we can plot this cropped data. We will show the AOI boundary on the same plot for scale.

ggplot() + 
  geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
  geom_raster(data = CHM_TUD_manual_cropped_df,
              aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  coord_sf()

Extract Raster Pixels Values Using Vector Polygons

Often we want to extract values from a raster layer for particular locations - for example, plot locations that we are sampling on the ground. We can extract all pixel values within 20m of our x,y point of interest. These can then be summarized into some value of interest (e.g. mean, maximum, total).

Image Source: National Ecological Observatory Network (NEON) To do this in R, we use the extract() function. The extract() function requires:

The raster that we wish to extract values from, The vector layer containing the polygons that we wish to use as a boundary or boundaries, we can tell it to store the output values in a data frame using df = TRUE. (This is optional, the default is to return a list, NOT a data frame.) . We will begin by extracting all canopy height pixel values located within our aoi_boundary_HARV polygon which surrounds the tower located at the NEON Harvard Forest field site.

tree_height <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), df = TRUE)

str(tree_height)
'data.frame':   621642 obs. of  2 variables:
 $ ID   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ layer: num  5.57 5.22 5.18 4.77 2.88 ...

When we use the extract() function, R extracts the value for each pixel located within the boundary of the polygon being used to perform the extraction - in this case the aoi_boundary_HARV object (a single polygon). Here, the function extracted values from 18,450 pixels.

We can create a histogram of tree height values within the boundary to better understand the structure or height distribution of trees at our site. We will use the column layer from our data frame as our x values, as this column represents the tree heights for each pixel.

ggplot() + 
  geom_histogram(data = tree_height, aes(x = layer)) +
  ggtitle("Histogram of CHM Height Values (m)") +
  xlab("Tree Height") + 
  ylab("Frequency of Pixels")

We can also use the summary() function to view descriptive statistics including min, max, and mean height values. These values help us better understand vegetation at our field site.

summary(tree_height$layer)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.460   0.000   0.375   4.343   8.523  36.729 

Summarize Extracted Raster Values

We often want to extract summary values from a raster. We can tell R the type of summary statistic we are interested in using the fun = argument. Let’s extract a mean height value for our AOI. Because we are extracting only a single number, we will not use the df = TRUE argument.

mean_tree_building_height_AOI <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), fun = mean)

mean_tree_building_height_AOI
         [,1]
[1,] 4.342554

It appears that the mean height value, extracted from our LiDAR data derived canopy height model is 22.43 meters.

Extract Data using x,y Locations

We can also extract pixel values from a raster by defining a buffer or area surrounding individual point locations using the extract() function. To do this we define the summary argument (fun = mean) and the buffer distance (buffer = 20) which represents the radius of a circular region around each point. By default, the units of the buffer are the same units as the data’s CRS. All pixels that are touched by the buffer region are included in the extract.

Image Source:National Ecological Observatory Network (NEON)

Let’s put this into practice by figuring out the mean tree height in the 20m around the tower location (point_HARV). Because we are extracting only a single number, we will not use the df = TRUE argument.

point_Delft <- st_read(here("data", "delft-leisure.shp"))
Reading layer `delft-leisure' from data source 
  `/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD New
mean_tree_height_tower <- extract(x = CHM_TUD,
                                  y = point_Delft,
                                  buffer = 20,
                                  fun = mean)

mean_tree_height_tower
  [1]          NA  0.18382418          NA          NA          NA          NA
  [7]  0.14741693          NA          NA  0.95534103          NA          NA
 [13]          NA          NA          NA          NA          NA          NA
 [19]          NA  1.08576819  1.77810457  1.50813700          NA          NA
 [25]          NA          NA          NA          NA  3.41709123          NA
 [31]          NA  5.44526759          NA          NA          NA          NA
 [37]          NA          NA          NA          NA          NA          NA
 [43]          NA          NA          NA          NA  1.26241786          NA
 [49]          NA          NA          NA          NA          NA  3.63058446
 [55]          NA          NA          NA  2.64284311          NA          NA
 [61]          NA          NA  2.18177335          NA          NA          NA
 [67]          NA          NA          NA          NA  1.71895883          NA
 [73]  4.59002702          NA  9.22261513  3.71037319  3.67939371  2.67243673
 [79]  4.01660438          NA          NA          NA          NA          NA
 [85]          NA          NA          NA  5.09733638          NA          NA
 [91]  3.60859259          NA          NA          NA  1.93923010  5.21877630
 [97]          NA          NA          NA          NA          NA          NA
[103]          NA          NA          NA  1.77163454          NA          NA
[109]          NA          NA  1.70136037  2.58267291  4.21806988          NA
[115]          NA          NA          NA          NA          NA          NA
[121]          NA          NA 10.30380493  2.28110616          NA          NA
[127]          NA          NA          NA  3.64850905          NA  0.08129702
[133]          NA  2.06486905 11.40954040          NA          NA          NA
[139]  2.08243861  1.11961589  6.40306065  6.41666083  6.71693856  5.15843022
[145]  4.27773571          NA          NA          NA          NA          NA
[151]          NA          NA          NA          NA          NA          NA
[157]          NA  3.32952196          NA          NA          NA          NA
[163]          NA          NA          NA          NA          NA          NA
[169]          NA          NA          NA          NA  2.07421889          NA
[175]  1.94442299  2.60783294          NA          NA          NA          NA
[181]          NA  7.78255215          NA          NA          NA  1.21956164
[187]  2.31018698          NA          NA          NA          NA  7.29981885
[193]  2.72441063          NA          NA          NA  0.04019622  6.68842712
[199]  6.17018350  1.59578616  0.66848060  5.40900358          NA          NA
[205]          NA  2.25253693          NA          NA  3.20768940          NA
[211]          NA  0.25329162          NA          NA          NA  1.00615888
[217]  3.88431955          NA 10.92242381          NA          NA          NA
[223]          NA          NA  1.77282882          NA          NA          NA
[229]          NA          NA          NA          NA  1.34524428  1.75811156
[235]  1.87248210  1.37764249  1.50010280  2.11123471  2.40555998  1.10978271
[241]  0.85010793  2.65174185  2.63599304  2.40464816          NA          NA
[247]          NA          NA          NA 10.16022778          NA          NA
[253]          NA          NA          NA          NA          NA          NA
[259]          NA          NA  0.61521977          NA  2.39372841          NA
[265]          NA          NA          NA          NA          NA          NA
[271]          NA          NA          NA          NA          NA          NA
[277]          NA          NA  5.96494518  2.29402613          NA          NA
[283]  2.52794784  0.73166368          NA  1.58280219  4.89034569  0.07626293
[289]  0.40960883  3.79256709          NA          NA          NA  3.07805207
[295]          NA          NA          NA          NA

Challenge: Extract Raster Height Values For Plot Locations

Solution

leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) %>% 
  filter(leisure %in% c("playground", "picnic_table"))
Reading layer `delft-leisure' from data source 
  `/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD New
# extract data at each plot location
mean_tree_height_plots_TUD <- extract(x = CHM_TUD,
                                       y = leisure_locations_selection,
                                       buffer = 20,
                                       fun = mean,
                                       df = TRUE)

# view data
mean_tree_height_plots_TUD
     ID       layer
1     1          NA
2     2  0.95534103
3     3          NA
4     4          NA
5     5          NA
6     6          NA
7     7          NA
8     8          NA
9     9          NA
10   10          NA
11   11  1.08576819
12   12  1.77810457
13   13          NA
14   14          NA
15   15          NA
16   16          NA
17   17          NA
18   18  3.41709123
19   19          NA
20   20          NA
21   21  5.44526759
22   22          NA
23   23          NA
24   24          NA
25   25          NA
26   26          NA
27   27          NA
28   28          NA
29   29          NA
30   30          NA
31   31          NA
32   32          NA
33   33          NA
34   34          NA
35   35          NA
36   36          NA
37   37  3.63058446
38   38          NA
39   39          NA
40   40          NA
41   41  2.64284311
42   42          NA
43   43          NA
44   44          NA
45   45  2.18177335
46   46          NA
47   47          NA
48   48          NA
49   49          NA
50   50          NA
51   51          NA
52   52  1.71895883
53   53  3.67939371
54   54          NA
55   55          NA
56   56          NA
57   57  1.93923010
58   58  5.21877630
59   59          NA
60   60          NA
61   61          NA
62   62          NA
63   63          NA
64   64  1.77163454
65   65          NA
66   66          NA
67   67          NA
68   68  1.70136037
69   69  4.21806988
70   70          NA
71   71          NA
72   72          NA
73   73          NA
74   74          NA
75   75          NA
76   76          NA
77   77          NA
78   78  2.28110616
79   79          NA
80   80          NA
81   81  3.64850905
82   82          NA
83   83  0.08129702
84   84          NA
85   85  2.06486905
86   86 11.40954040
87   87          NA
88   88          NA
89   89          NA
90   90  2.08243861
91   91  1.11961589
92   92  6.40306065
93   93  6.41666083
94   94  6.71693856
95   95          NA
96   96          NA
97   97          NA
98   98          NA
99   99          NA
100 100          NA
101 101          NA
102 102          NA
103 103  3.32952196
104 104          NA
105 105          NA
106 106          NA
107 107          NA
108 108          NA
109 109          NA
110 110          NA
111 111          NA
112 112          NA
113 113          NA
114 114          NA
115 115          NA
116 116          NA
117 117          NA
118 118          NA
119 119  1.94442299
120 120          NA
121 121          NA
122 122          NA
123 123  7.78255215
124 124          NA
125 125          NA
126 126  1.21956164
127 127          NA
128 128          NA
129 129  7.29981885
130 130  2.72441063
131 131          NA
132 132  6.17018350
133 133  5.40900358
134 134          NA
135 135          NA
136 136          NA
137 137          NA
138 138          NA
139 139          NA
140 140  3.88431955
141 141          NA
142 142          NA
143 143          NA
144 144          NA
145 145          NA
146 146          NA
147 147          NA
148 148          NA
149 149  1.34524428
150 150  1.75811156
151 151  1.87248210
152 152  1.37764249
153 153  1.50010280
154 154  2.11123471
155 155  2.40555998
156 156  1.10978271
157 157  0.85010793
158 158  2.65174185
159 159  2.63599304
160 160  2.40464816
161 161          NA
162 162          NA
163 163 10.16022778
164 164          NA
165 165          NA
166 166          NA
167 167          NA
168 168          NA
169 169          NA
170 170          NA
171 171          NA
172 172  2.39372841
173 173          NA
174 174          NA
175 175          NA
176 176          NA
177 177          NA
178 178          NA
179 179          NA
180 180          NA
181 181          NA
182 182          NA
183 183          NA
184 184          NA
185 185          NA
186 186          NA
187 187          NA
188 188  5.96494518
189 189  2.29402613
190 190          NA
191 191  2.52794784
192 192  0.73166368
193 193          NA
194 194  3.07805207
195 195          NA
196 196          NA
197 197          NA
198 198          NA
# plot data
ggplot(data = mean_tree_height_plots_TUD, aes(ID, layer)) + 
  geom_col() + 
  ggtitle("Mean Tree Height at each Plot") + 
  xlab("Plot ID") + 
  ylab("Tree Height (m)")


Summary and keypoints.

We have seen how to crop a raster to the extent of a vector layer and how to extract values from a raster that correspond to a vector file overlay. In short: - Use the crop() function to crop a raster object. - Use the extract() function to extract pixels from a raster object that fall within a particular extent boundary. - Use the extent() function to define an extent.